clip-02-05.jl
Load Julia packages (libraries) needed for the snippets in chapter 0
using StatisticalRethinking, Optim
#gr(size=(600,600));snippet 3.2
p_grid = range(0, step=0.001, stop=1)
prior = ones(length(p_grid))
likelihood = [pdf(Binomial(9, p), 6) for p in p_grid]
posterior = likelihood .* prior
posterior = posterior / sum(posterior)
samples = sample(p_grid, Weights(posterior), length(p_grid));
samples[1:5]5-element Array{Float64,1}:
0.727
0.489
0.607
0.68
0.721snippet 3.3
Draw 10000 samples from this posterior distribution
N = 10000
samples = sample(p_grid, Weights(posterior), N);10000-element Array{Float64,1}:
0.859
0.567
0.417
0.701
0.738
0.562
0.73
0.719
0.618
0.513
⋮
0.596
0.594
0.59
0.652
0.684
0.759
0.436
0.319
0.551In StatisticalRethinkingJulia samples will always be stored in an MCMCChains.Chains object.
chn = MCMCChains.Chains(reshape(samples, N, 1, 1), ["toss"]);Object of type Chains, with data of type 10000×1×1 Array{Float64,3}
Iterations = 1:10000
Thinning interval = 1
Chains = 1
Samples per chain = 10000
parameters = toss
2-element Array{ChainDataFrame,1}
Summary Statistics
. Omitted printing of 1 columns
│ Row │ parameters │ mean │ std │ naive_se │ mcse │ ess │
│ │ Symbol │ Float64 │ Float64 │ Float64 │ Float64 │ Any │
├─────┼────────────┼──────────┼──────────┼────────────┼────────────┼─────────┤
│ 1 │ toss │ 0.636631 │ 0.137942 │ 0.00137942 │ 0.00131646 │ 10405.3 │
Quantiles
│ Row │ parameters │ 2.5% │ 25.0% │ 50.0% │ 75.0% │ 97.5% │
│ │ Symbol │ Float64 │ Float64 │ Float64 │ Float64 │ Float64 │
├─────┼────────────┼──────────┼─────────┼─────────┼─────────┼─────────┤
│ 1 │ toss │ 0.349975 │ 0.545 │ 0.645 │ 0.737 │ 0.878 │
Describe the chain
MCMCChains.describe(chn)2-element Array{ChainDataFrame,1}
Summary Statistics
. Omitted printing of 1 columns
│ Row │ parameters │ mean │ std │ naive_se │ mcse │ ess │
│ │ Symbol │ Float64 │ Float64 │ Float64 │ Float64 │ Any │
├─────┼────────────┼──────────┼──────────┼────────────┼────────────┼─────────┤
│ 1 │ toss │ 0.636631 │ 0.137942 │ 0.00137942 │ 0.00131646 │ 10405.3 │
Quantiles
│ Row │ parameters │ 2.5% │ 25.0% │ 50.0% │ 75.0% │ 97.5% │
│ │ Symbol │ Float64 │ Float64 │ Float64 │ Float64 │ Float64 │
├─────┼────────────┼──────────┼─────────┼─────────┼─────────┼─────────┤
│ 1 │ toss │ 0.349975 │ 0.545 │ 0.645 │ 0.737 │ 0.878 │
Plot the chain
plot(chn)
snippet 3.4
Create a vector to hold the plots so we can later combine them
p = Vector{Plots.Plot{Plots.GRBackend}}(undef, 2)
p[1] = scatter(1:N, samples, markersize = 2, ylim=(0.0, 1.3), lab="Draws")
snippet 3.5
Analytical calculation
w = 6
n = 9
x = 0:0.01:1
p[2] = density(samples, ylim=(0.0, 5.0), lab="Sample density")
p[2] = plot!( x, pdf.(Beta( w+1 , n-w+1 ) , x ), lab="Conjugate solution")
Add quadratic approximation
plot(p..., layout=(1, 2))
End of 03/clip-02-05.jl
This page was generated using Literate.jl.